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ABSTRACT 

A description is given of the algorithms implemented in the AstroBEAR adap- 
tive mesh refinement code for ideal magnetohydrodynamics. The code provides 
several high resolution, shock capturing schemes which are constructed to main- 
tain conserved quantities of the flow in a finite volume sense. Divergence free 
magnetic field topologies are maintained to machine precision by collating the 
components of the magnetic field on a cell-interface staggered grid and utilizing 
the constrained transport approach for integrating the induction equations. The 
maintenance of magnetic field topologies on adaptive grids is achieved using pro- 
longation and restriction operators which preserve the divergence and curl of the 
magnetic field across co-located grids of different resolution. The robustness and 
correctness of the code is demonstrated by comparing the numerical solution of 
various tests with analytical solutions or previously published numerical solutions 
obtained by other codes. 
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Introduction 



The development of efficient and accurate numerical algorithms for astrophysical flow 
has become of great interest to the astrophysical community. Variable resolution approaches 
have provided an avenue for efficient simulation of hydro dynamical flow including multi- 
physical effects which involve substantial variation in length scale. Adaptive Mesh Refine- 
ment (AMR) has been recognized as one of the most versatile and efficient approaches to 
enable the simulation of multi-scale phenomena for which fixed-resolution simulation is either 
impractical or impossible. AMR discretizations employ a hierarchy of grids at different levels. 
High resolution is applied only to those regions of the flow which would otherwise be subject 
to unacceptably large truncation error. The utility of the AMR approach is underscored by 
the extensive list of codes that are tar geted toward astrop hysical r esearch which utiliz e AMR. 



he list includes AstroBEAR, En zo fO'She a et a.. 



(Truelove et al 



( IFromang et al. 



i998; Klein 1999; Crockett et a 



2006^ . RIEMANN jBalsaralboOlh and AMRVAC fceppens et al. 



2004T) . Flash flFrvxell et alJl2000h . Orion 
2005h . Nirvana Jzierierl 12005a). Ra mses 



20031 ) and 



the list of codes for whic h the develo pment AMR capabili ty is in progress including Athena 
rtGardiner fc Stone! l2005h and Pluto dMignone et~al1(2007h . 



The simulation of magnetized flow is of particular interest to astrophysical researchers 
owing to the utility of numerical magnetohydrodynamics (MHD) in modeling a wide range of 
astrophysical phenomena. The leading line of recent research in this area has focused on the 
appli cation of higher order Godunov methods to numerical MHD flRyu &: Joneslll995l ; lBalsara 
19981 ). The conservative formulation and proper upwinding employed by these methods 
enable accurate simulation of strongly supersonic flow. Because of this unique capability, 
such methods are often referred to "high resolution shock capturing" (HRSC) methods. 

While HRSC methods have long been recognized as the defacto standard for the simula- 
tion of supersonic hydrodynamical phenomena, their popularity among researchers interested 
in magnetized flow has been slowed because standard HRSC approaches to MHD fail to main- 
tain the solenoidality constraint on the magnetic field (V • B = 0). If not corrected, local 
divergences in the magnetic field arising from this short coming usually grow rapidly, causing 
anomalous magnetic forces an d unphysical plasma trans port which eventually destroys the 
correct dynamics of the flow fjBrackbill fe Barnes! Il980l ). Early practitioners of numerical 
MH D therefore relied h eavily on finite difference methods as employed by codes such as 
Zeus IStone &: Norman! (119921 ) which maintain t he solenoidality constraint exactly despite 
their inferior shock capturing ability (iFalld 120021 ) . Later works focused on improved HRSC 
which either eliminate the development of divergences in the magnetic field, or mitigate the 
effect of local divergence errors on the dynamics of the flow. In one such approach, a pro- 
jection operator is devised, usually by solving a Poisson equation, which removes numerical 
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divergences from the grid after each time- step (jBalsaralll998l ; I Jiang &: Wulll999l ; iKim et al. 
19991 ; IZachary et al.lll994t iRyu et al.lll995l ). The primary li mitation of this a pproach is that 
non-trivial boundary value problems become i ndeterminate ( Rvu et al.lll998T) . In the second 
so called "8-wave" approach first explored by iPowell et al.l (119991 ) . alternative formulations 
of the MHD equations are constructed to prevent the local build up of magnetic divergence 
by advecting monopoles to other re gions of the grid wher e they are of less consequence to the 
dynamics of the flow. The work of iDedner et al.l (120021 ) augments this approach by adding 
source terms to the system which act to counter the effect of local divergence in the mag- 
netic field on the dynamics of the flow. A third approach known as constrained transport 
utilizes a multidimensional, divergence-preserving update procedure for the magnetic field 
compone nts which are collated on a staggered grid centered on the compu t ational volume in- 



terfaces. (lEvans fc Hawlev 



1988 



the most accurate results in the tests of 



Rvu et al.lll99a Balsara fc Spicerl Il999l : iDai fc Woodward 



1998; Londrillo & Del Zanna bood : Izieglerll2005bh^ This a pproach has been show n to provide 



Tothl fboooh and balsara fc Kiml tooi ). 



The combination of AMR spatial discretizations with HRSC would seem a natural choice 
in order to satisfy the desire for a high accurate, computationally efficient and versatile strat- 
egy for the simulation of magnetized plasma flow. The implementation of the aforementioned 
adaptations to HRSC methods for MHD in an AMR framework however, poses several chal- 
lenges. Divergence cleaning schemes utilizing a Poisson projection operator are ill-suited for 
AMR applications owing to difficulty in handling the p rojection step along internal bound- 



aries on a patchwork of grids at different resolutions. IPowell et al.l ( 19991 ) found that the 
application of their 8-wave method on an AMR grid hierarchy, local divergence errors on one 
level of the AMR hierarchy caused local divergence of comparable magnitude on all levels. 
The main drawback of this method in an AMR context is that the unphysical effects of 
local diverg ences in the magnetic field are not diminished by the application of additional 
refinement. ICrockett et al.l (120051 ). however, have constructed an approach suitable for AMR 
applications which combines an approximate projecti on operator w i th the dive rgence advec 



tion a nd dampening the effects of local divergence of IPowell et al.l (119991 ) and IDedner et al. 
(l2002h . 



Retaining the divergence-free property of the solution obtained through the applica- 
tion of the constrained transport approach on hierarchical grids requires application of a 
divergence-preserving prolongation operator which interpolates the magnetic field from a 
coarse mesh a fine mesh, a divergence-preserving restriction operator which maps the fine 
mesh magnetic field onto a coarser mesh and furthermore requires that the evolution of 
the magnetic field be consistent between coll ocated meshes of different resolution. Two ap- 
proaches to these challenges ha ve emerged. iBalsaral (l200lf) has generalized the divergence 
free reconstruction procedure of iBalsara fc Spicer ll999h to devise a prolongation operator 
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based on a piece-wise qu adratic i n terpo lation procedure that is divergence preserving in the 



RIEMANN MHD code. iLi fc Lil (120041 ) present an adaptation of Balsara's pr ocedure that 



simpli fies its implementation for problems involving arbitrary refinement ratios. iToth &: Roe 



( 120021 ) have devised a prolongation operator by solving an algebraic system which enforces 
the maintenance of the volume average curl and divergence between grids of different reso- 
lution in the AMR hierarchy. 

In this paper we provide a concise description of the algorithms and tests of the As- 
troBEAR HRSC AMR MHD code. AstroBEAR is comprised of several numerical solvers, 
integration schemes, and radiative cooling modules for astrophysical fluids. The code's AMR 
capability is derived from the AMR engine of the BEARCLAW boundary embedded adaptive 
mesh refinement package for conservation laws. This code utilizes the constrained transport 
approach to adapting HRSC methods to the MHD system of equations. To our knowledge , 



AstroBEAR is the first AMR code to utilize the prolongation operator of IToth fc Rod (120021 ) 
to maintain the V • B = constraint. By combining multi-physics capabilities relevant to 
simulation astrophysical plasma flow, AMR, and a wide selection of HRSC integration proce- 
dures, AstroBEAR will serve as a valuable research tool. The authors intend that this paper 
will serve as a reference for future works that apply the code and provide a concise recipe for 
robust, reliable and accurate HRSC solution strategies for MHD on AMR grid hierarchies. In 
§2]we describe the several HRSC schemes and divergence preservation strategies implemented 
in the code. In $3] we provide an overview of the AMR algorithm, highlighting the stages 
which require special treatment of the magnetic field. In §Hwe provide a concise description 
of the prolongation, restriction and coarse to fine refluxing procedures required to preserve 
the divergence and consistency of the magnetic field across an AMR hierarchy of grids. In 
§[5] we comment on the results of several test and example problems with particular emphasis 
on the relative strengths and weaknesses of the various HRSC schemes implemented in the 
code. In £0we provide a synopsis and discussion of the main results of the paper. 



2. Numerical Method 

AstroBEAR provides an adaptive mesh refinement framework for the integration of 
conservation laws of the form, 

|Q + ^F.(Q) + ^F y (Q) + ^F,(Q) = S(Q). (1) 
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In this work we focus on the equations of ideal MHD which are written in the conservative 
form as: 
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with gas density p, velocity v, total volumetric energy density £, thermal pressure P, mag- 
netic field B, and electric field E. In the proceeding system of equations, we have chosen 
units for the electric and magnetic field in the plasma so that factors of 4n do not appear in 
the equations. The last three equations in the system follow from Faraday's law, 

^B + VxE = 0, (3) 

implicit to which is the constraint that initially solenoidal magnetic field topologies remain 
solenoidal, 

V-B = 0. 

The equations are brought to a closed form via Ohm's law for a perfectly conducting medium, 

E = -v x B (4) 
and the polytropic equation of state for an ideal gas, 

P=( 7 -l)(£-pv 2 /2-B 2 /2). (5) 



In the remainder of this section, we describe the details of the shock capturing numerical 
schemes available in our code to integrate the solution to equations of the form of equation 
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[TJ and modifications thereof which ensure the solenoidal constraint on the magnetic field. 
The purpose of this description is to provide a concise and complete illustration of the steps 
necessary to build the code. The details of any of the particular solution strategies may 
be obtained by consulting the original work credited for the particular strategy. For a more 
pedagogically oriented re yiew of h i gh res oluti on shock cap turing schemes, we refer the reader 
to the excellent books of iLevequd (120021 ) and iTorol (119991 ). In the remainder of this section, 
we consider only the solution method for the homogeneous part of the conservation law with 
S = 0. The effect of non-zero source terms which may be used to includ e additional physics 



are ha ndled using the operator splitting technique described in §2.1 of ICunningham et al. 



(120051 ). The micro-physical source terms included in the code to model the effects of radia- 
tive cooling, and ti me dependent, non-equilibr ium ionization and H2 chemistry are cataloged 
in the appendix to ICunningham et al.l (120051 ) . In appendix A we catalog the MHD source 
terms and modifications to the numerical scheme that are employed in applications involving 
cylindrical axisymmetric flow. We adopt the notation that superscripts denote time, the first 
subscript denotes the direction of vector components, and the last three subscripts denote 
the spatial location on the computational grid. The superscript is omitted from temporally 
varying quantities which are to be evaluated at time t. Equations which demonstrate oper- 
ations that take the same form in each of the x, y, and z directions are written only for the 
x-direction sweep. In these cases, the y sweep can can be recovered by replacing y — > x and 
j — > i and the z sweep can can be retrieved by replacing z — > x and k — *■ i. 

Numerical integration of the system of conservation laws is achieved using the finite 
volume method. The basis of the finite volume quadrature can be realized by discretizing 
the integral form of equation [TJ yielding the following unsplit procedure for advancing the 
conserved field forward in time by an increment At: 
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where F", F™, and F" are suitably accurate numerical approximations to the inter-cell flux, 
spatially averaged over the inter-cell area and temporally averaged from t to t+At. Construc- 
tion of higher order schemes in more than one dimension can, in some cases, be simplified 
by utilizing a direction-split approach 
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QjJ.fc ~~ + ^ (-^x3,i-l/2j,fc(Q!j',fe) ^ 1 x3,i+l/2J,fe(QiJ,fc)J ■ CO 

In equation ([7]) we have used the superscript in parenthesis to denote intermediate states 
between direction sweeps. The ordering of the component directions (xi,X2,xs) cycles over 
different permutations of the three coordinate directions on successive time-steps. In two 
dimensions the ordering is [{x,y), (y,x)] and in three dimensions we have found the most 
consistent results by cycling over both the cyclic and anti-cyclic permutations: 

[(x, y, z), (z, x, y), (y, z, x), (z, y, x), (x, z, y), (y, x, z)] 

The subject of the following subsections is the procedure to compute the numerical flux 
which is comprised of three steps 1) reconstruction (interpolation) to zone edges §2.11 2) 
upwinding the solution of the Riemann problem at each zone edge §2.2} and 3) temporal 
evolution of the field of conserved quantities §2.31 Our code implements several methods for 
carrying out each of these steps which may be utilized in the combination that best tailors 
the integration strategy to the requirements of the application at hand. 



2.1. Spatial Reconstruction 

We construct a spatially second order accurate integration procedure via suitable re- 
construction of the state in each computational cell from the volume average state within 
that cell and its neighbors. We define a "primitive variable" operator, L p (Qj) = Pj = 
[p, v x , v y , v z , E — pv 2 /2 — B 2 /2, B x , B y , B z \ , which converts the conserved state variables 
into a form more suitable for interpolation. Interpolation of the primitive variables, rather 
than the conserved, has the advantage that the reconstructed state at grid edges are guar- 
anteed to have non-negative pressure. We write the reconstruction from the volume-average 
state variables collated at grid centers to the grid interface as: 

Pl,i+1/2 = Pj + 2^+'* 

PjM-i/2 = P * - \<P~,i- (8) 

Figure [1] shows a cartoon schematic of the volume average field of primitive variables on a 
one dimensional grid (solid lines), the cell reconstruction (dotted lines) and the location of 
the left and right restricted states. 
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\v2 X i + 1 



Fig. 1. — Schematic of diagram of the volume average field of primitive variables on a one 
dimensional grid (solid lines), the cell reconstruction (dotted lines) and the location of the 
left and right restricted states, Pl and Pr. 
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The grid interface conservative fields are then computed as 



Ql,?:-i/2 

Q-R,i-l/2 



: L p (Pi,i_i/ 2 ) 

Lp 1 (PRJ-l/2) 



(9) 



The code implements a user selectable choice of three different seco nd order s patia l 
reconstruction methods. The first is the MUSCL reconstruction method of I Van Leer! (119791 ) 
using a slope limiter to maintain monotonicity, 



fr ±>i = LIMITER (Pi - P f _!, P m - P 4 ) 



(10) 



Many limiter functions may be constructed. AstroBEAR implements a choice of three lim- 
iters, which operate on the state vector parameters in equation [TD] in a component-wise 
fashion. In order of decreasing levels of diffusion introduced into the scheme, these limiters 
are the "min-mod" limiter, 



MINMOD(x,y) 
the limiter of vanLeer 
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MIN(ABS(x), ABS(j/))SIGN(x) otherwise 
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o x v 

x+y 



if xy < 
otherwise 



(12) 



and the "monotonized-centered" limiter 





MC(x,y) 



(13) 



if xy < 

MIN(2x, \(x + y),2y) otherwise 

The second reconstruction met hod is the l o cal hy perbolic harmonic variation of the piecewise 
hyperbolic method (PHM) of IMarquinal (11994J ). The PHM reconstruction prescribes the 
components of (f>±i as 



if \5 L \ < 10" 14 and \6 R \ < 10" 14 

d Ax r]± otherwise 



(14) 
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h 
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if \h\ < 10 
if \5 R \ < 10- 14 
otherwise, 



or (SlSr > and £r < 5^) 
or (5 L 5 R > and S L < S R ) 



P - Pi- 



Ax 



P 



i+1 



Ax 



Note that r\ has a removable singularity about k = with lim K ^ +(^) = lim K ^ - {if) = 1. 
When evaluated with 8-byte precision, the expression for rj begins to diverge from its true 
solution for \k\ < 10~ 5 and we set r\ — > 1 in this region using the piecewise expression given 
above. 

The third method is the reconstruction procedure of the piecewise parabo lic method 

(PPM ) of lColellal dl984h. A detailed overview of the PPM method is availab le from lMiller k Colella 
( 120021 ). and iMignone et al.l (120051 ). Appendix B of iMignone et al.l ( 120051 ) already provides 
a concise description of the PPM reconstruction procedure which we do not repeat here. 
We note that in the numerical examples using the PPM reconstruction pre sented later in 



this p aper, we have taken a different approach to numerical oscillations than IMignone et al. 
( 120051 ). In particular, we maintain monotonicity via the MINMOD limiter (equation [TT 



rather than the more compressive Van Leer limiter and we do not included the dissipation 
mechanisms of §B.l of IMignone et al.l (120051 ). 



2.2. Upwinded Numerical Flux 



The numerical flux is computed by upwinding the waves associated with the Riemann 
problem defined by Ql,«-i/2 on the left and Q_R,i_i/2 on the right of each computational cell 
interface. AstroBEAR implements three differe nt methods for computing the upwinded flux, 
the HLLD flux as described in as described in lMiyoshi fc Kusanol (120051 ) the Roe flux, and 
the Marquina flux. 



The Roe flux method iRod (119861 ) calls for the decomposition of the cell edge states into 
eigenmodes of a linearized sy stem matrix. I n the case of MHD we utilize the approximate 
linearized Riemann solver of iRyu fc Jonesl (119951 ) . We write this decomposition in terms 
of the eigenvalues a m ( Q ) , right eigenvectors R m (Q) and left eigenvectors L m (Q) given in 
§2.2 of IRyu fc Jonesl (119951 ) where the subscript m denotes the m th eigenmode. Singularities 
which arise, in certain limi ting cases, in the norm alization of the eigensystem are avoided in 
the manner prescribed by iRoe fc Balsaral (119961 ) . The inter-cell flux across cell boundaries 
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at constant x is computed as 
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(15) 



m=l 



where the eigendecomposition is carried out about a suitable average of the states to the left 
and right of the cell interface which we denote as < Ql,i-i/2, Qi?,i-i/2 >• 
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The flux functions of iDonat fc Marquinal (119961 ). first proposed by lShu fc Osher Jl989h 
for scalar equations, take the form 

f?x,i-l/2 = El=i «ill/2 R m,i-l/2 + a £-l/2 R m,i-l/2 ( 17 ) 
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ft 



i-1/2 



k |ELiMAx(|a^_ 

r o 



i i 

1/2 I' l a m,i-l/2 



l) L m,i-i/ 2 Qi,i-i/2 otherwise 



if a m,i-l/2 a £,i-l/2 > 

and °m,i-l/2 > 
if a ™,i-l/2 a m,i-l/2 > Rl8) 
aIld <i-l/2 < 



a 



R 

i-1/2 



El=l L m, ! -l/2 F (Qfl,'-l/2) 

k |ELiMAx(l< 4 



1/2 1' \ Uj m,i-l/2 



if a m,i-l/2 G m, i -l/2 > 

and a L m ._ ll2 > 

if ffl i,i-l/2 (l m, ! -l/2 > ((19) 
aIld <i-l/2 < 



Lf,i_i/ 2 QL,i-i/2 otherwise 



The numerical flux across cell boundaries at constant y, F y j_i/ 2 , and z, F Zj fc_i/2, are com- 
puted in the analogous manner. In this case a separate eight decomposition of the system 
matrix on each side of the cell interface is required. 
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(20) 
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In the coarse of carrying out simulations in the hydrodynamic limit, we have found the 
Marquina flux to be advantageous for modeling certain astrophysical phenomena (discussed 
below). However, the method fails to converge for some MHD shock tube test problems that 
give rise to compound structures, that i s, structures com posed of a shock and rarefaction of 
the same wave family moving together (jBrio fc Wulll988l). In particular , the method fails to 
converge for the Riemann problem discussed in §5 of iRyu fc Joned (119951 ) . We have overcome 
this deficiency by constructing an adaptation of Marquina's flux which introduces a small 
a moun t of nu merical diffusion by utilizing the same averaging procedure as the flux formula 
ofQ Jl986h . 
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(21) 



For the methods calling for the decomposition of a linearized approximation to the 
system matrix (the Roe flux and o ur ad aptation of the Marquina flux), we use the arithmetic 
average method of IRyu fc Joned (119951 ) where the system is linearized at the cell interface 

as 



< P 



RJ 



-1/2 >— I- 1 L, i— 1/2 A R, i-1/2) 

Numerical tests have shown that averaging the net of the magnetic and thermal pressure 
rather than the thermal pressure alone provides more accurate results. Therefore we set 



f-pRJ 
\^L, i-1/2 



)RJ 

R, i 



i)/2- 



(22) 



tRJ 



[p, v x , v y , v z , P + B 2 /2, B x , B y , B z ] 



(23) 



In the puerly hydrodynamic limit the code uses the average state of iRod (119811 ). This 
averaging guarantees desirable property that 



(9F(< Q L|i -i/ 2 , Qrj-i/2 >) 

dQ(Q L ~ Qr) 



F,(Q L )-F x (Q /? ) 



(24) 



which is not true, in general, for the arithmetic ave r age li nearization. We therefore employ 
the arithmetic average linearization of IRyu fc Joned (119951 ) for MHD applications and revert 
to the Roe average linearization only for purely hydrodynamic applications. 

In general, the Roe flux provides the least diffusive fo r mulat ion. This is because the flux 
formulations based on the method of iDonat fc Marquinal (119961 ) revert to the more diffusive 
local Lax-Freidrichs upwinding for transonic eigenmodes. This additional component of 
numerical diffusion is is advantageous for simulating certain astrophysical phenomena. It 
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is sufficient to pre vent the developme nt of rarefaction shocks without the need to introduce 
an "entropy fix" (IHarten et al.lll976l ). It is also sufficient to dam pen the development o f 
carbuncles, even-odd decoupling and related numerical pathologies ([Sutherland et al.ll2003l ) 
that are particularly problematic in grid-aligned flows which exhibit strong radiative cooling. 



2.3. Temporal Reconstruction 

The temporal update operation for the grid centered values is given by equation [6j 
Replacing n — > t and utilizing any of the methods for computing the numerical flux in 
the previous section achieves an integration that is first order in time. We implement four 
different methods to obtain second order temporal accuracy by performing this update using 
time-centered estimates of the numerical flux. 

The MUSCL-Hancock predictor-corrector temporal discretization achieves second order 
accuracy by advancing the grid-face interpolated states by a half-time-step using a one 
dimensional predictor step. The predictor step is carried forward according to: 

Qn,t% = Qk^i/2 + ^(F,(Qk,-i)-F,(Qi, m )) + ^S(Q„). (25) 

Note that the predictor step uses the cell centered, volume average fluxes. The Riemann 
problem at the cell interfaces is not solved and the upwinded flux at the cell-faces are not 
needed for this step. 

The MUSCL-Hancock corrector step calls for the construction of a second order, time 
centered numerical flux which is computed by applying any of the upwinding procedures of 
( §2.2p using Q^^irt, and as ^ ne anc ^ right states for the Riemann problem at 

the i — 1/2 cell interface. We note the update to the time centered state described here, use 
of this procedure in more than one dimension requires the application of the operator split 
integrator to retain the second order accuracy of the scheme for multidimensional flow. The 
fully second order accurate update is therefore carried out via application of equation [7] with 
the time centered flux with n=dt/2 

The application of the predictor-corrector schemes like the MUSCL-Hancock approach 
described above will in some cases necessitate the application of a protection procedure to 
ensure pressure and density positivity of the predictor interface states: 

Pl^% - max(p^ 1 / / 2 2 ,io- 2 min(p1 i1 _ 1/2 + p^_ 1/2 ),io- 14 ) 

- MAX(p^/; 2 ,10- 2 MM^ 
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Plti/l - MAX(^ A ; / / 2 2) 10- 4 MIN(Pi il _ 1/2 + P^-i/ 2 ) ) 10- 14 ) 
Pr£X% <~ MAX(P^ A f; 2 ,10- 4 MIN(P^_ 1/2 + P^i/ 2 ),10- 14 ). (26) 

The second temporal integratio n opt i on im plemented in the code is the two step Runge- 
Kutta temporal update operator of IShul (119881 ). In the first step F* i _ 1 ^ 2 and Q*~^_i/ 2 are 



computed via the application of a first order update step given by equation[6]with n = At and 
any operator split microphysical effects (e.g., source terms). In the second step F(Q*"t_* , 2 ) 
is computed by a second application of the spatial interpolation and upwinding procedure 
on the grid of Qjj,"^* / 2 data. The second order, time centered fluxes are then computed via 
the interpolation formula, 

= \ (F(QU/ 2 ) + F(QS A ! /2 )) • (27) 

This flux may be used to carry forward a fully second order, unsplit update via equation 
O The unsplit nature of the Runge-Kutta time stepping is a significant advantage of the 
method. The Runge-Kutta scheme also has a comparative advantage in that it is pressure 
positivity-preserving and therefore more robust. However, the method entails somewhat 
greater computational cost due to the solution of twice the number of Riemann problems 
per grid cell per time step as the MUSCL-Hancock approach. In addition the scheme suffers 
a somewhat restrictive time-step stability condition. The maximum numerically stable time- 
step At is computed in terms of the Courant condition as 

At < MAX[ C m ^ hk ) + MAx C m ' l V l/2 ' k ) + MAx C m,t f~ l/2 ) for all i,j, k\. (28) 
LV Ax ' K Ay ' Az 

In practice, we estimate the next time-step increment from the maximum wave speed en- 
countered during the preceding integration sweeps as: 

i± /~i n t a tr a irr/^ m ^- 1/2,7, fc a m,i,j — l/2,k ^m,i,j,k— 1/2 \ £ n • ■ ; 1 / nn \ 

dt next = CFL MAX[( — , — " J , — '—) for all k). 29 

Ax Ay Az 

where CFL is a user tunable parameter. We typically choose, CFL ~ 0.8 for one dimensional 
calculations and CFL ~ 0.4 for multidimensional problems. Future revisions of the code 
will inc l ude th e multidimensional corner transport upwind (CTU) reconstruction method of 
Colellal Jl990h and enhan cements to this method for the integration of the MHD equations by 



Gardiner &: Stone! (120051 ). By explicitly including the effect of transverse-propagating waves 
at each grid interface, the CTU scheme retains numerical stability for larger time-steps, 
CFL < 1, while capturing greater accuracy. 
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2.4. Constrained Transport Scheme 

While Faraday's law (equation [3]) guarantees solenoidality of dB/dt, and therefore the 
maintenance of solenoidal magnetic field topologies, the Gudonov-based conservative field 
update procedures described in the previous section do not provide any such guarantee. This 
is because no provision has been made in the construction of the integration procedure that 
would enforce the divergence constraint on the magnetic field. Each of the fluxes used in the 
conservative update procedure, are second order approximations to the exact area averaged 
fluxes where 



F x + 0(Ax 3 
F y + 0(Ay 3 
F z + Q(A^ 3 



(30) 



Therefore, the divergence of the magnetic field after a conservative update will also contain 
high order truncation errors, with VB = 0(Ax 3 )+0(Ay 3 ) + 0(Az 3 ). Local departures from 
V • B = usually grow rapidly, causing anomalous magnetic forces an d unphysical plasma 
trans port which eventually destroys the correct dynamics of the flow (jBrackbill fc Barnes 
19801 ). Two strategies have emerged for adapting Godunov-based MHD schemes so that the 
divergence-free constraint is exactly maintained. In the first approach a projection operator 
is devised , usually by s o lving a Poisson equation which removes numerical divergences from 
the g rid faalsaral fl~998l : Ijiang fc Wulll999l : iKim et ali[l~99~9l : IZacharv et all 11994 : Ir^u et al. 



19951 ). Solving the Poisson equation is somewhat computationally expensive and particu- 
larly algorithmically complex on AMR grid hierarchies. The second approach utilizes a more 
multidimensional approach toward the numerical quadrature of Faraday's law by utilizing a 
conservative formulation of Stoke's theorem to represent magnetic field components at stag - 
gered collocation points (IBalsara fc Spicerl Il999l ; iDai &: Woodward Il998l; iRyu et al.lll998l ). 
Following the nomenclature first commissioned by Evans fc Hawley (jl988 ). this approach is 
commonly referred to as "constrained transport" (CT) in the literature. The AstroBEAR 
code utilizes the CT approach to maintain a divergence free fiel d, primarily due to the limi- 
tations of the first approach in AMR applications. Furthermore, IBalsara &: Kiml (120041 ) have 
demonstrated the superiority of the staggered grid approach in the context of a stringent as- 
trophysically motivated test problem involving the interplay of strong shocks with radiative 
cooling. 

The basis of the constrained transport approach is realized by applying of Stoke's the- 
orem to Faraday's law and integrating over each face of a control volume. This yields an 
expression for the face-average normal component of dB/dt at each control volume interface. 
Spatial discretization of the line integral around the control volume interfaces calls for the 
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average electric field parallel to each edge of the control volume. Temporal discretization 
reveals an explicit update procedure for the normal component of the magnetic field at each 
computational volume interface. The resulting discretized equations, 

r>t+At ni , ^ 

D x,i-l/2,j,k - n x,i-l/2,j,k "I" AyAz 

(^yEy t i-i/2j,k+i/2 — ^yEy t i-i/ 2 ,j,k-i/2— 

AzE Z) i-i/2 t j+l/2,k + Az£' 2 , i _i/2j-l/2,it) 
r>t+At jyt i ^£ 

^W-iA* - n vM-W + AxAz 

(AzE Z) i+l/2 ) j-l/2,k — AzE Z) i-\/2,j-l/2,k~ 

AxE Xt ij- 1 / 2t k+l/2 + AxE Xti j^i/2 t k-l/2) 

r>t+At ni , ^ 

n z,i,j,k-l/2 ~ n z,i,j,k-l/2 "I" AxA2/ 

[AxExjj+xfak-i/z — AxE x ^j- 1 /2,k-i/2— 
AyEy ii+1 / 2 j t k-i/2 + AyEy^i/ 2 j t k-i/2) ■ (31) 

gives the desired CT update operator. Note that the centered difference discretization of the 
divergence of the magnetic field, 

/•n d\ Bx,i+l/2,j,k ~ Bx,i-l/2,j,k , By,i,3+l/2,k ~ ^y,i,j-\l2,k . -D«,i,i,fc+l/2 — -Dz,iJ,fc-l/2 

(V ' B) ^ = Ax + Ay + Az 

(32) 

is preserved 

(V • B)£* = (V • B)^ (33) 

and the CT update operator maintains a divergence free representation of the magnetic field 
provided that the initial magnetic field is divergence-free. We emphasize that the CT update 
procedure requires that the components of the magnetic field be collated at the center of the 
zone faces to which they are orthogonal and that the component of the electric field parallel 
to each computational cell interface be known. The spatial location of the desired electric 
and magnetic field components are illustrated in figure [2J 




Fig. 2. — The location of the staggered electric and magnetic field components on a com- 
putational cell centered at position k. 
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Specification of the CT update procedure is completed via a suitable construction of 
the required electric field components at each grid edge. The fluxes of the MHD equation 
(equation [2]) can be expressed in terms of the electric field using Ohm's law (equation @J. 
The numerical inter-cell fluxes computed during the conservative update step described in 
the preceding subsections, provide a second order, shock capturing approximation to the 
components of the electric field at the center of each computational cell interface. The CT 
update scheme, however, calls for the electric field at the grid edges. A reconstruction of the 
cell-face electric fields to the cell edges that retains second order accuracy is given by: 



E 



tt+At/2 7 / /.c-t-iAr/z ^c^-iAr/z /.c-^-iAr/z ri^Lii/A 

x,i,j-l/2,k-l/2 ~ W z 7,i,j,k-l/2 7,i,j-l,fe-l/2 _ Jy 8,i,j-l/2,k ~ Jy 8,i,j-l/2,fc-l 



(t+At/2 



et+At/2 



t+At/2 



pt+At/2 



E 



t+At/2 

y,i-l/2,j,k-l/2 



, ( nt+At/2 nt+At/2 nt+At/2 nt+At/2 

8,i-l/2,j,k ' fx 8,*-l/2,j,fc-l ~ J z 6,i,j,k-l/2 ~ J z 6,i-l,j,fc-l/2 



rpt+At/2 
2,i-l/2,j-l/2,fc 



Wy 6,t,i-l/2,Ai ' Jy 6,i-lJ-l/2,fc J-x 7,i-l/2,j,k J-x 7,i-l/2,j-l,fc I ' \ 6q: ) 



P t+At/2 



t+At/2 



ft+At/2 



where the cell-face electric fields have been writ ten in terms of ce ll inter face fluxes. Setting 
/* = i 71 * and k — 1/4 recovers the CT scheme of iBalsara fc Spicerl (Il999l ). Many procedures 
for averaging the cell interfa ce flux component s to construct the cell-edge electric field com- 
ponents. The CT scheme of iRyu et al.1 (119981 ) can be expressed in the form of equation [M] 
by retaining only the advective part of the inter-cell flux, 



fx 6,7,8 = V x [0, B y , B Z ] T 

% 6,7,8 = v y [B x , 0, B Z ] T 
f 2 6,7,8 = v z [B x ,B y ,0] T 



(35) 



and setting k — 1/2. The upwinded, time centered, cell-face fluxes ft +At ^ 2 are constructed 



from the cell centered flux components /* according to the same procedure given in §2.21 
This update procedure retains the upwinding abilities of the conservative update scheme 
while maintaining a divergence free solution to the magnetic field via the correspondence 
of the components of the magnetic flux with the components of the electric field called for 
by the CT update procedure. Construction of the electric field from the components of 
the intercell numerical flux therefore has the advantage of int r oduci ng very little numerical 
dissipation into the solution. Of the schemes tested by iTothl (120001 ) . those that utilize this 
approach to reconstructing the electric field at grid edges produce the most accurate results. 
One implementation detail that should be noted is that the update of the components of the 
magnetic field collated along the boundary of the computational domain requires the flux 
components that are parallel to that boundary extending into the first row of ghost cells 
along the boundary. The conservative update procedure must therefore extend one row of 
computational cells into the ghost cell region during integration sweeps transverse to the 
boundary even though no conservative update is applied to the boundary cells. Extension 
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into the ghost region ensures that the components of the numerical flux required in the CT 
update step are computed. 

At the end of each CT update step, the volume centered cell solution for the magnetic 
field, computed during the Gudonov update, is discarded in favor of the solution provided 
by the CT update. At this stage in the algorithm, the cell centered the magnetic field is 
recomputed according to the procedure: 

B x ,i,j,k — ^ {Bx,i+l/2,j,k + B xi _i/2j y k) 
1 

By,i,j,k = ~ (-By,i,i+l/2,fe + ^2/,i,i-l/2,fc) 

B z ,i,j,k — 2 (B z ,i,j,k+l/2 + ^z,i,j,k-X/2) ■ (36) 

Even though the solution of the magnetic field is advanced in time at cell interfaces, we retain 
the grid of volume average magnetic field components in order to compute the magnetic 
pressure at the cell centers. Also, the volume average magnetic field and total energy are 
synchronized with the cell-face magnetic field after every CT update step in order to preserve 
the volume averaged thermal energy. 

R 2 R 2 

°i,j,k < °i,j,k 2 2 

The total energy is thereby adjusted such that the CT update preserves thermal energy. This 
optional step avoids numerical difficulty with negative thermal pressure that arise in strongly 
magnetized problems at the expense of energy conservation. We view this as an acceptable 
trade-off, particularly for astrophysical applications that are not energy conserving due to 
radiative energy losses. 

Because the normal component of the magnetic field is known at cell faces, it is not 
necessary to interpolate these values during the spatial reconstruction step of §2.11 Recalling 
that we have defined the 6th, 7th and 8th components of the primitive field vector as the x, 
y and z components of the magnetic field (equation [8]), we set 
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(38) 

in lieu of the spatial reconstruction procedure when using an unsplit scheme (equation [H]). 
With a direction split scheme (equation [7]) the partial update of the cell-centered field from 
prior direction sweeps is interpolated back to the cell edge and added to the cell edge state 



-20- 



during the spatial reconstruction step (§sr) in order to retain second order accuracy with 
multidimensional flow 
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(39) 



where the superscript (a — 1) refers to the state after the previous direction sweep. While this 
procedure does maintain the second order accuracy achieved via direction splitting it does 
suffer from the flaw that the reconstructed states are not exactly divergence free. Instead, 
the reconstructed states given by MUSCL-Hancock time stepping procedure described in 
§2.31 and any direction split scheme suffers from the undesirable property that, 



-1/2; 
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is not maintained exactly. Also note that the CT update procedure is applied after each 
stage of the multistage Runge-Kutta temporal update procedure such described in §2.31 



We generally find the divergence error introduced via the split reconstruction procedure 
is small and that because we apply a divergence free CT update for the magnetic field, the 
effect of this small error does not grow rapidly. The effect of this error appears as inexact 
ev olution of the z-compone nt of the magnetic field in the magne t ic fiel d loop advection test 
of iGardiner fc Stond (120051 1 as shown in §5.21 iGardiner Stond (120051 ) in troduced a two di- 
mensi onal unsplit scheme and later a three dimensional extension thereof IGardiner fc Stone 
( 120081 ) where they devise an unsplit reconstruction that is not subject to t his error at the 
expense of considerably increased algorithmic and computational complexity. iBalsaral (120041 ) 
has shown that unsplit MHD schemes are less diffusive in some tests than the dimensionally 
split counterpart. Because of this we recommend the dimensionally unsplit Runge-Kutta 
scheme which does provide for exactly divergence free reconstructed states for AtroBEAR 
MHD applications. Never the less the direction sp lit MUSCL-Hancock upda te has proven to 
be re l iable in earlier pu rely hydro dynamic works ( Cunningham et al. 2006al lbl; iDennis et al. 



20081 ; lYirak et al.l 120081 ) . and for this reason we leave this option available for MHD appli- 
cations. We note that the MUSCL-Hancock scheme is efficient in terms of computational 
cost, requiring only one Riemann solver per grid cell per ti me step while ret aining stability 
for CFL < 1 and that exhaustive testing of this and other (IRyu et al.l 1 19981 ) direction split 
MHD schemes have shown good results. 
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2.5. Summary or Numerical Methods 



The menu of options available in our code includes three spatial reconstruction methods 
§2.11 linear (MUSCL), piecewise hyperbolic (PHM) and piecewise parabolic (PPM), two 
temporal reconstruction methods §2.31 : MUSCL-Hancock and Runge-Kutta, four different 



Riemann solvers / upwinding procedures §2.21 the HLLD flux, the Roe flux, the Marquina 
flux, and an adaptation of the Marquina flux that is better suited for magnetized flow involv- 
ing compound wave structures and two different constrained transpo rt schemes for preserving 
the solenoidality co nstraint on magn etized flows §2.41 the method of iBalsara fc Spicerl (119991 ) 
and the method of iRyu et al.l (119951 ). The user of the code may choose the method for each 
of these operations which are summarized in table [TJ The advantage of this approach is that 
a given simulation may be performed using several different methods. The solution strategy 
which is optimal for the physical regime of a given simulation may be readily applied. 
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Table 1: Numerical Method Options. 



Spatial Reconstruction 


Temporal Reconstruction 


Flux Function 


CT Scheme 


MUSCL 


MUSCL-Hancock 


Roe 


Ryu et al. 


PHM 


Runge-Kutta 


Marquina 


Balsara & Spicer 


PPH 




Adapted Marquina 








HLLD 
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3. Adaptive Mesh Refinement 

The central feature of the BEARCLAW framework on which AstroBEAR is based is 
that it provides a framework for Adaptive Mesh refinement (AMR). Under AMR, regions of 
the flow that are susceptible to large discretization errors are carried forward on a computa- 
tional grid of higher resolution while flow features not requiring high resolution for adequate 
numerical convergence are carried forward on a computational grid of lower resolution. Two 



app roaches to AMR have em erged: (1) the block-based method of iBerger fc Oligerl (119841 ) 



and IBerger fc Colellal (119891 ) which constructs a patchwork of refined grids that optimally 
covers all of the cells on the next-coarsest level that are heuristically identified for refinement 
(2) an al ternative approa ch where individual computational cells are refined or derefined sep- 



arately (iKhokhlovl Il998l ) . Our code employs the former approach. In this section we give 
an overview of the AMR algorithm and stages of the AMR algorithm that require special 
attention in handling grid face magnetic field components. For discussion we use the term 
"parent grid" to designate an underlying block on the next coarser level, "parent level" to 
designate all of the grid blocks that are one level coarser, "child grids" to designate those grid 
blocks that are one level finer and "child level" to refer to grids that are of higher resolution 
by one refinement ratio than the current grid. Advancement of an AMR hierarchy of grids 
is carried out according to the pseudo-code algorithm given in appendix [Bl A schematic of 
the update procedure is shown in figure EH for an AMR hierarchy of three levels of refinement 
where each level has a refinement ratio of 2. Curved horizontal arrows represent integration 
of all grids on a given refinement level. The algorithm is adaptive in time, with each level 
advanced in time increments dti eve i = dt\ eve i^\jr where r is the refinement ratio of the level. 
Gray vertical arrows represent restriction and refluxing of the solution on refined grids to 
their parent level. Black vertical arrows represent prolongation of the solution from the 
coarse level to its parent level. 



-24- 




Fig. 3. — A schematic of the update procedure for an AMR hierarchy of three levels of 
refinement where each level has a refinement ratio of 2. Curved horizontal arrows represent 
integration of all grids on a given refinement level. Gray vertical arrows represent restriction 
of refined grids to their parent level. Black vertical arrows represent prolongation of the 
solution from the coarse level to its parent level. 
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The first step of the update procedure, "Set Ghost," calls for the initialization of the 
ghost cells that are exterior to all grids on the given level. Figure H] shows an example of 
an AMR hierarchy containing one grid on the root level (level=0), and one refined level 
(level=l). The interior of each grid is delineated by solid lines and the extended grid which 
include the interior and ghost region of each grid are delineated by dashed lines. The face- 
centered magnetic field components that coincide with boundaries delineating the interior of 
each grid are treated as interior cells. We classify ghost cells into three categories: 1) same- 
level ghost cells that coincide with the interior of another grid on the same level appear white 
in the figure, 2) physical ghost cells that lie outside of the computational domain appear 
dark gray in the figure, 3) child-level ghost cells that coincide with the interior of grids on 
the parent level appear light gray in the figure. Same-level ghost zones are initialized to 
the state of the interior of the coincident grid. Physical ghost cells are initialized according 
to user specified boundary conditions, the code provides three physical boundary options; 
constant extrapolation, reflecting or periodic. The initialization of parent-level ghost cells is 
carried out as a part of the "Grid Adapt" procedure which will be discussed below. 
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Fig. 4. — An example of an AMR hierarchy containing one grid on the root level (level=0), 
and one refined level (level=l). The interior of each grid is delineated by solid lines and 
the extended grid which include the interior and ghost region of each grid are delineated by 
dashed lines. 
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The "Grid Adapt" procedure determines the arrangement of and initialization of a new 
AMR hierarchy on a given level that tracks the evolution of flow features in the solution. 
A user-specified truncation error estimation procedure is applied to each grid that is one 
level coarser. AstroBEAR uses either the maximum absolute value of the primitive vector 
gradient or Richardson extrapolation (IBerger fc Qligerlll984l ; iBerger fc Colellal 119891 ) as the 
options for error estimation. Zones on the parent grids with estimated error greater than 
a user-specifi ed tolerance are flagged fo r refinement. We employ the patch- wise clustering 
algorithm of IBerger fc Rigoutsosl (119911 ) to determine arrangement of refined grid patches 
that optimally overlays all of the zones flagged for refinement. For easier parallel implemen- 
tation, we require that each grid has only one parent. In figure |5l the new arrangement 
of grid patches on the first refinement level with interior boundaries and ghost boundaries 
delineated by solid and dashed lines respectively. The region coinciding with the previous 
arrangement of AMR patches are shaded in gray. Interior regions of the new grid arrange- 
ment that coincide with interior regions of the previous patchwork of grids on the same level 
are initialized by copying the field values from the previous grids. The previous patchwork 
of grids on this level are then released from memory. The ghost zones, and interior zones 
that do not coincide with the interior of the previous grid patches are initialized from the 
parent grid via a prolongation operator. The code prolongs the cell-centered conserved fields 
using interpolation from the parent grid. The cell-face grid of magnetic field components 
are initialized using a divergence-preserving prolongation operator in order to preserve the 
integrity of the CT update procedure. We will discuss this operator in detail in §4.21 
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Fig. 5. — The new arrangement of grid patches on the first refinement level with interior 
boundaries and ghost boundaries delineated by solid and dashed lines respectively. The 
region coinciding with the previous arrangement of AMR patches are shaded in gray. 
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The "integrate" step advances the solution on each grid on the given level using one of 
the integration procedures discussed in $2] as specified by the user. The time-adaptive nature 
of the AMR engine imposes the difficulty that boundary information from parent grids is not 
available during each step of child grid integration cycle. As shown in figure[3]for a refinement 
ratio of two, child ghost zones are temporally synchronized with their parent grids only every- 
other time cycle. To accommodate this, we incorporate an "extended" arrangement of ghost 
cells. Each refined grid carries a strip of ghost zones that extends the width r x m& c cells 
beyond the interior of the grid where r is the refinement ratio, and rribc is the number of 
ghost cells required by the integration stencil. The MUSCL-Hancock reconstruction operator 
( §2.ip without direction splitting requires two ghost cells, yielding rribc = 2. An extra ghost 
zone is required in multidimensional problems when using direction splitting and constrained 
transport in order to compute the flux components necessary to compute the EMF at the 
grid corners, yielding m& c = 3. The Runge-Kutta temporal integration method ( §2.31) has 
been implemented via additional rows of ghost cells to fully update the interior region with 
rribc = 4. The update of the extended region of ghost cells is illustrated in figure [6] for 
a refinement ratio r = 2. The initial representation of the field is at time t. The first 
integration cycle carries all cells interior to the first ring of rribc ghost zones shaded in dark 
gray forward in time from t to t + dt/2. On the second integration cycle, the interior cells 
shaded in dark gray are carried forward in time from t to t + dt/2. The second integration 
step can be carried out because the outermost ring of rribc ghost zones are outside of the 
domain of influence on the interior cells. 




Fig. 6. — A schematic of the arrangement of the extended ghost region on a grid with a 
refinement ratio of 2. The initial representation of the field is at time t. The first integration 
cycle, illustrated in the left panel, carries all cells interior to the outer ring of ghost 
zones shaded in dark gray forward in time from t to t + dt/2. On the second integration 
cycle, illustrated in the right panel the interior cells shaded in dark gray are carried forward 
in time from t to t + dt/2. 
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The "Synchronize to Parent" call carries forward the synchronization of the solution on 
a given level to its parent grid. The synchronization procedure is composed of two steps: 
application of coarse to fine level refluxing and restriction of the solution. In the refluxing 
step, a spatio-temporal average of the flux from child grid interfaces along coarse / fine grid 
boundaries is compared to the flux as computed during the integration of the coarse grid. 
A correction is applied to the parent grid so that the effective flux across this boundary is 
equal to that as computed on the child grid. The restriction step calls for the coarsening 
and injection of field data from the interior of the child grids into their respective parent 
grids. We employ volume weighted average restriction of cell centered conserved fields. The 
restriction of face-centered magnetic field components requires special attention to maintain 
the divergence-free property of the magnetic field which will be discussed in the next section. 



4. Divergence Preserving Restriction and Prolongation Operators 



In §2.41 we discussed the importance of satisfying the V • B = in order to maintain the 
integrity of the numerical solution to the MHD equations and demonstrated an adaptation 
to Goduonov-based schemes that preserves the divergence of the magnetic field throughout 
the calculation. The use of AMR imposes an additional challenge that must be overcome 
to maintain the divergence of the magnetic field throughout the calculation. Specifically, 
the restriction step, which maintains the consistency on coarse levels with the finer levels of 
refinement in the grid, and the prolongation step, which initializes coarse representations of 
the solution to finer grids in regions that have been flagged for refinement must not introduce 
divergence errors into the solution. The volume average and bilinear interpolation procedures 
utilized for the restriction and prolongation of the cell-centered conserved fields cannot be 
adapted to operate on the grid-edge magnetic fields in a divergence pres erving m a nner. Two 
approaches to the divergence preserving prolongation h ave emerged. iBalsaral (20011 ) has 



generalized the divergence free reconstruction procedure of iBalsara &: Spicer 



(1999|) to devise 



a piece-wise quadratic interpolant that is divergence preserving. iLi &: Lil (120041 ) present an 
adaptation of Balsara's proc edure that simplifies its implementation for problems involving 
arbitrary refinement ratios. iToth fc Rod (120021 ) has devised a prolongation procedure by 
solving an algebraic system which enforces the maintenance of the volume average curl and 
divergence over a coarse grid cell. 

In this subsections that follow, we present the restriction and prolongat ion formulae 
empl oyed in our code. In particular, we present the prolongation procedure of (IToth &: Roe 
20021 ) in a less compact form than that of the original authors which are more readily 
transcribed into computer code. To simplify implementation our code only allows refinement 
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ratios of two between successive levels of refinement. 



4.1. Restriction 



Figure H illustrates the locations of the magnetic field components utilized by the con- 
strained transport update procedure where a refined grid (denoted with dotted lines) has 
been embedded within one coarse grid cell (denoted with solid lines) for a two dimensional 
calculation. We denote the magnetic field components on coarse grid faces as (3, on refined 
cell faces that are coincident with coarse grid faces (the exterior faces) as b, and on refined 
cell faces that do not coincide with coarse grid faces (the interior faces) as B. We have 
adapted a short hand notation where the first, second, and third character in the superscript 
to the refined grid magnetic field components either represents the location above (+), below 
(-) or aligned with (0) the center of the coarse cell in each of the x, y, and z directions. 
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x, i 1/2, j 



b". . 




b + + 



X, i+1/2, j 



Ay 



Fig. 7. — Illustration of the locations on a two dimensional grid of the staggered electric and 
magnetic field components with a refined grid, delineated with dotted lines, that has been 
embedded within one coarse grid cell, delineated by solid lines. 



We require that the fine-to-coarse grid synchronization step maintain the cell interface 
magnetic field in an area-averaged sense, 





5 (^ + + bt) 



(41) 



in two dimensions, and 



Py,i,j±l/2,k 



Px,i±l/2,j,k 



Pz,i,j,k±l/2 



\(bt ±+ + by ±+ + b^ + b y ^) 
~ (bt +± + K +± + 6^ ± + K~±) 



(42) 



in three dimensions. Equations UU & H21 must be satisfied while preserving the divergence 
of the magnetic field along coarse/fine grid boundaries. Simultaneous satisfaction of these 
properties along coarse/fine grid interfaces and this restriction condition is accomplished by 
applying a suitable restriction of the electric field from fine to coarse grids before performing 
the CT update (equations EJJ) on the coarse level. Manipulation of equations EH subject 
to the constraint provided by equations HJJ & I421 at all times yield the desired electric field 
restriction operator as 



in two dimensions and 

pt+dt/2 l i t+dt/4 t+3dt/4 t+dt/4 t+Zdt/4 

2,M-l/2,fc-l/2 ~~ 4V e x,«-l/4j-l/2,fc-l/2 e a;,i-l/4,j-l/2,fe-l/2 e x,i+l/4,j-l/2,A;-l/2 e x,j+l/4,j-l/2,fc-l/2 

rpt+dt/2 i / t+dt/4 t+3dt/4 t+dt/4 t+Mt/4 

^y,i-l/2,j,k~l/2 ~ ~4~\ e y,i-l/2,j-l/4,k-l/2 ' e j/,i-l/2,j-l/4,fc-l/2 ' e y,j-l/2,j+l/4,fc-l/2 + e y,i-l/2,j+l/4,k-l/2 

pt+dt/2 _ i / t+dt/4 t+Zdt/4 t+dt/4 t+Zdt/4 a 

z,i-l/2,j-l/2,k ~ 4V e z,i-l/2J-l/2,A;-l/4 e z,i-l/2j-l/2,fc-l/4 e 2 ,i-l/2J-l/2,fc+l/4 e z,i-l/2,j-l/2,fe+l/4v 

in three dimensions where E is the electric field on the coarse grid and e is electric field of 
the refined grid. In figure [3] note that level is advanced by an increment from t to t + dt/2 
using the electric field computed at time t + dt/2 in one step while the next child level, level 
1, is integrated by the same increment in two steps using the electric field at time t + dt/4 
to integrate from t to t + dt/2 and electric field at time t + 3dt/4 to integrate from t + dt/2 
to t + dt. The temporal averaging of the electric field is necessary due to the temporal 
refinement capability of the code. The time averaging in equations H3] & HH ensures that 
evolution of the magnetic field across level boundaries remains divergence free and consistent 
in the sense of equations HU & 1421 



E 



,t+dt/2 

'«,t-l/2J-l/2 



1 ( t+dt/4 
2\ e -z,i-l/2,j-l/2 



+ e 



t+J.dt/4 x 
: 2 ,i-l/2,i-l/2J' 



(43) 
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4.2. Prolongation 

The prolongation step initializes a newly refined grid from its parent grid that is coarser 
by one level of refinement. The prolongation of the face centered magnetic field is carried 
out in two steps. In the first stage, the exterior faces that coincide with the edges of an 
already refined grid block are set by copying the field values from the coincident face of the 
already refined grid. The exterior faces that do not coincide with the edges of an already 
refined region are computed via a bilinear interpolation of the coarse representation of the 
field given by 

b x + = Px,i±l/2,j + -jbyfix 
bx~ = Px,i±l/2,j - ^yflx 
= A/,jj'±l/2 + -j^xPy 
h y ± = Py,i,j±l/2-^ x f3 y , (45) 

in two dimensions and 



b± ++ 


= Px,i±l/2,j,k + 2< 


SyPx + SzPx) 


b± + - 


= Px,i±l/2,j,k + 


SyPx - SzPx) 


bt + 


= Px,i±l/2,j,k + 


Syfix + SzPx) 


bt- 


= Px,i±l/2,j,k + ^ ( 


SyPx - SzPx) 


b + y ±+ 


= Py,i,j±l/2,k + 2< 


SxPy + S z (3 y ) 


^ +± - 


1 

= (3y,i,j±l/2,k + ^( 


8 x (3y - 8 z (3y) 


by ±+ 


= Py,i,j±l/2,k + ^( 


-S x y - 5 z (3 y ) 


bt~ 


= Py,i,j±l/2,k + ^( 


-5 x [3 y - 5 z (3 y ) 


b ++± 

z 


= Pz,i,j,k±l/2 + ~( 


8 X Pz + 5y(3 z ) 


b +-± 

"z 


= Pz,i,j,k±l/2 + ^{ 


8 X Pz - Sy(3 z ) 


K +± 


= Pz,i,j,k±l/2 + -( 


-5 x f3 z + 5 y (3 z ) 




= Pz,i,j,k±l/2 + 


-8x(3 z ~ SyPz), 
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in three dimensions, where we compute the spatial jumps using the same slope limiter (equa- 
tions [TT]{T3]) as the base scheme: 



&x@* — LIMITER(/3* j j + ij_i/2,ifc-l/2 ~~ P*,i,j-l/2,k-l/2, /3*,i,j-l/2,k-l/2 ~ P*,i-lj-l/2,k-l/V 
8yP* = LIMITER(/3* j i_l/2,i+l,fc-l/2 — fi*,i-l/2,j,k-l/2, (3*,i-X/2,j,k-\/2 — (3*,i-X/2,j-l,k-l/2) 
8zP* = LIMITER(^* j j_i/2j'-l/2,Jfe+l — fl*,i-l/2,j-l/2,k, P*,i-l/2,j-l/2,k — P*,i-l/2,j-l/2,k-0'7) 

In the second stage of the prolongation procedure we interpolate from the exterior faces 
of the refined grid cell to the interior faces. In three dimensions, the 12 refined inter-cell faces 
that are interior to the coarse cell are constructed via an interpolation from the magnetic 
field components collated at t he 24 refined grid e xterior face centers that coincide with a 



coarse cell interface. Following IToth fc Rod (120021 ) we construct this interpolation to satisfy 



the constraint that it be divergence and curl preserving. For readability and conciseness, 
we illustrate the steps required to derive such an interpolation procedure in detail only for 
the two dimensional case. We then present the analogous solution in three dimensions. 
The volume average divergence, computed via the application of Gauss's law around the 
perimeter of the exterior faces is: 

d 00 = l -{bt + + 5+" - K + - b-- + &++ + b y + - b+- - b~) (48) 

where we have used the bar accent to denote division by the discretization width in the 
normal direction on the coarse level (see figure Ej), i.e. b x = b x /Ax, and b y = b y /Ay. 
The curl of the field is constructed by bilinear interpolation of the requisite exterior field 
components to the origin of the refined cell as: 



c 



00 



A V n++ T-+ , T+- 7 — \ Ax 



(fe++ - b y + + 6+" - b-~) - —(bt + - bf + b~ + - b--). (49) 



2 Ax v y y • y y > Ay 

The divergence centered at each of the four refined cell interiors is: 

D- = 2(B°- - b~~ + B-° - b; 



D+- =2(5+--S£- + S+°-5+-) 

ZT+ = 2(B° X + - b~+ + b y + - B y °) 

D ++ =2(6++-S£+ + 6+ + -B+°) (50) 



and the curl at the origin of the refinement cells implied by the interior field values is: 

The desired interpolation procedure is determined by imposing the condition that each of 
the refined grid cells contributes equally to the divergence of the refinement region, d 00 = 



C° z ° = - B y °) - ^{Bl + ~ B° x -). (51) 
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D = D~ + = D + ~ = D ++ , and that the curl implied by the interior faces equals the 
curl interpolated from the exterior faces, c° z ° = C°°. Because only three of four divergence 
conditions are linearly independent, we have a system of four independent linear equations 
for the four desired interior field values. The solution to the system, in matrix form, is: 



b ex t 
Bint 

[A] 



hit 



[K- 

2 

1 
1 

A] bext- 



5-0 
V 



2 

2 2 

-1 -1 1 

-1 -1 1 



1 -1 



1 

2 




2 






2 



i 
i 
o 

2 



(52) 



In three dimensions, we follow the same procedure to derive the solution for twelve 
interior refined cell face fields from 24 exterior refined cell face fields. In particular we 
have eight equations (seven linearly independent) for the volume integrated divergence over 
each of eight refined cell interiors, d 000 = D = D + = D~ + ~ = D~ ++ = = 
D + ~ + = D ++ ~ = D +++ , and six equations (five linearly independent) by computing the 
three components of the curl, each at two positions relative to the center of the coarse cell, 

c x ~ { ~ / x\-dx/4 — { ^x\+dx/i} c y — ^yl-dy/i ~ ^y\+dy/<ii dIlu c z ~ ( - y z\-dz/i ~ ( - y z\+dz/i- 

solution of this system for the desired interior field values yields the desired prolongation 
procedure which can be written in the form of equation [52] as: 

bext = [ b. 



Bint = 



— b— + b-+- &-++ b+— b + - + b ++ - b +++ 

3j X £C 3j X CC 3j 

~ b y -+ 6-++ 6+- b; s ' 6++- b+++ 

~ K- + K + - K ++ K- bt- + 5++- bf++] T 

b :°- b : o+ b -^ B t-° B t +o f 
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5. Numerical Tests, Examples & Discussion 



The AstroBEAR numerical sche mes have been tested against a suite of test problems, 
including; the one dimensional tests of iRyu fc Jonesl (119951 ) and the two dimensional tests of 
Ryu et al.l (119951 ). Except for the failure of the original Marquina flux formulation to con- 
verge for problems involving compound MHD wave structures as discussed in §2.21 each of 
the methods recover results that are equivalent to those of earlier researchers barring minor 
differences due to different levels of numerical truncation error. In the remainder of this 
section we present the results of some of these tests. The results of a few popular hydrody- 
namic shocktubes are presented which illustrate the differences between and limitations of 
each of the reconstruction and upwinding methods. In these cases we provide commentary 
which is intended to provide guidance as to the optimal choice of methods given the expected 
physical regime of a par ticular simulation . We also reproduce the results of several of the 
two dimensional tests of IRyu et al.l (119951 ) using AMR. These tests demonstrate the success 
of the divergence preserving restriction and prolongation operators presented in §H 



5.1. One Dimension 



Figure [8] shows the density re sulting from the ISodl (119781 ) shock tube problem denoted 
as "test problem 1" in the book by iTorol (119991 ) using several spatial reconstruction methods. 
The initial left and right Riemann problem states, number of computational zones and final 
time of the solution are presented in table [2J In all three cases the Runge-Kutta temporal 
integration and the flux formulation of Roe was used. We note that the result of this test 
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is not very sensitive to the choice of temporal reconstruction method or flux formulations. 
The shock tube consists from left to right of a backward propagating rarefaction, a contact 
discontinuity, and a forward propagating shock wave. The numerical diffusion caused by 
the truncation error of each of the spatial reconstruction methods is readily apparent in the 
figure. We characterize the diffusion of the method by the width of the contact discontinuity. 
The three methods are presented in order of increasing diffusion. The least diffusive method, 
PPM resolves the contact across three zones, and the most diffusive method, MUSCL resolves 
the contact across 8 zones. We note that the diffusion of the PPM and MUSCL methods 
may be improved by choosing a more compressive slope limiter than the MINMOD limiter 
used here at the expense of introducing oscillations near sharp discontinuities in the flow. 
Naturally, those methods which exhibit the least numerical diffusion are the most accurate. 
However, the PPM and PHM methods tend to "overshoot" the solution near sharp disconti- 
nuities, thereby introducing small oscillations into the solution. Such oscillations may drive 
numerical pathologies in simulations of flow with very low pressure near discontinuities such 
as persistently negative pressures. We note that some researchers have managed to reduce 
such oscillations by introducing addit ional sources of diffusion to the base scheme (see, for 
example §B1 of iMignone et al.l (120051 ) ). Simulation of astrophysical phenomena involving 
strong radiative cooling is particularly susceptible to this problem. In such cases the opti- 
mal solution strategy is determined as a trade-off between the desire for simultaneously low 
numerical diffusion and low oscillation. 
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Fig. 8. — Result of the ISodl (119781 ) shocktube (test problem 1 of iTorol (119991 )) using 
MUSCL/VL slope limiter (left), PHM (center) and PPM (right) spatial reconstruction. 
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Table 2: Sod 


Shocktube Parameters. 

left, x < 0.5 right, x > 0.5 


P 

P 

7 


1 0.125 



1 0.1 
1.4 1.4 


grid zones 
CFL 

t final 


128 

0.9 (0.5 for PPM) 
0.25 
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Figure [9] shows the result of the "1-2-0-3" strong expansion s hock tube of lEinfeldt et al. 



(Il99ll ). denoted as "test problem 2" in the book by iTorol (119991 ) . The initial left and right 
Riemann problem states, number of computational zones and final time of the solution is 
presented in table [3J This shock tube launches two rarefaction waves, one propagating to 
the left and one propagating to the right. The upper three panels show the test results 
utilizing the MUSCL - VL slope spatial and Ruge-Kutta temporal reconstruction methods 
with different flux upwinding procedures. This test problem was designed to illustrate the 
failure mode common to linearized Riemann solvers evident by the anomalous oscillation 
in velocity and anomalous rise in specific thermal energy (E t h/p) about the center of the 
grid in panel 9a. These anomalies are caused by the addition of small amounts of thermal 
energy to the solution in regions where the pressure becomes negative. We note that specific 
thermal energy anomalies such as these can result in a cascade of numeric al pathologies 



i n mu lti-physics simulations involving temperature dependent micro-physics. lEinfeldt et al. 



(119911 ) showed that the more diffusive HLLE Riemann solver, unde r certain conditions , guar - 
antees pressure positivity and therefore reduces such anomalies. iGardiner fc Stond (120051 ) 
have also demonstrated success applying the HLLE solver only in those regions where the 
linearized solver produces negative pressure or density. In panel 9b we show that application 
of the Marquina flux formulation also resolves the anomalous behavior. In §2.2l we presented 
an adaptation of the Marquina flux formulation which is better suited to magnetized flow 
involving compound wave structures. The result of this calculation utilizing the adapted 
flux formulation (panel 9c) shows that the adaptation retains the desirable features of the 
Marquina flux for this test problem. 
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x 



Fig. 9. — R esults of the lEinfeldt et al.l (I199ll ) "1-2-0-3" strong rarefaction shocktube (test 
problem 2 of iTorol (119991 )) using the Roe flux, the Marquina flux, and the adapted Marquina 
flux from top to bottom. 
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iauie o. Ji/in 


'e\dt Shocktube Parameters. 




left, x < 0.5 


right, x > 0.5 


P 


1 


1 




-2 


2 


P 


0.4 


0.4 


7 


1.4 


1.4 


grid zones 


128 




CFL 


0.8 




t final 


0.15 
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In figure [9] we reproduce the MHD shock tube test from figure 2a of iRyu fc Jones 

( 119951 ). Following the notation of IRyu fc Jones! (119951 ) . we denote the orientation angle of 
the magnetic field, \1/ = ta.n~ 1 (B z / B y ). The initial left and right Riemann problem states, 
number of computational zones and final time of the solution is presented in table [3J This 
Riemann problem demonstrates the ability of the code to correctly capture a number of 
MHD discontinuities. The discontinuities in the flow from left to right are, (1) fast shock, 
(2) rotational discontinuity, (3) slow shock, (4) contact discontinuity, (5) slow shock, (6) 
rotational discontinuity, and (7) fast shock. We note that the PHM spatial reconstruction 
e mployed for th i s test generally results in a level of truncation error slightly lower than that 
of IRyu &: Jonea (119951 ) with a slightly more diffusive result about of rotational discontinuities 
due to th e lack of a rotationa l discontinuity steepening procedure like that of equations 2.96 
- 2.98 of IRyu fc Joned (119951 ). We note that when applied to this particular problem, the 
MUSCL and PPM spatial reconstruction methods result in lower and higher truncation error 
respectively, in a manner that is similar to Sod problem discussed earlier. 




Fig. 10. — MHD shocktube result using PHM spatial reconstruction, MUSCL-Hancock tem- 
poral reconstruction, and the Roe flux function. 
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MHD Shocktube Parameters. 





left, x < 0.5 


right, x > 0.5 


p 


1.08 


1 




V x 


1.2 







Vy 


0.01 







V z 


0.5 







P 


0.95 


1 




B x 




2/v 


/ 4tt 


By 


3.6/V47T 


4/v 




B z 


2/VIn 


2/v 




7 


5/3 


5/3 




grid zones 
CFL 

t final 


512 
0.8 
0.2 
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5.2. Two Dimensions 

In this section we present the results of several two dimensional simulations to demon- 
strate the robustness of the divergence preserving scheme in AMR applications. The simu- 
lations presented in this section which utilize the AMR functionality of the code apply an 
AMR hierarchy of 3 levels using a refinement ratio of 2 between levels. 

We begin with the problem of the advectiq n of a weak magnetic field lo op following 



the prescription given in lGardiner fc Stone! (120051 ) and lGardiner fc Stond (120081 ). The initial 



conditions given in table [5] The initial face-centered magnetic field is generated from the 
analytic prescription of the vector potential (B = V x A) at cell corners using a centered 
difference formula on a grid extending from < L x < 2, < L y < 1 with resolution 
2N x N. In the left column figure [11] we show magnetic field lines (red) over a greyscale 
representation of the current density J = V x B for the initial condition (top) and after 
the advection of the loop through the periodic domain at t — 1 for both the Runge-Kutta 
(center) and MUSCL-Hancock temporal integration (bottom) with a resolution of N = 128. 
In both cases monotized-center limited linear spatial reconstruction with the flux function 
of Roe were used. The problem contains an initially singular current sheet along the surface 
and a singular spike at the center which is very sensitive to error in the evolution of the 
magnetic field. The propagation of these singular features through the grid serves as a very 
stringent test of the MHD update algorithm. We find that both methods maintain the 
correct location of these singularities and maintain magnetic field contours that are smooth 
and nearly symmetric about the current spike. 

The top and center portions of the right column of figure [TT] show the evolution of 
the spatially averaged magnetic energy normalized to the initial analytic value B for the 
MUSCL-Hancock and Runge-Kutta integration approaches respectively for different grid 
resolutions. We find that both methods show comparable accuracy and the expected con- 
vergence properties for this quantity. However, the direction split MUSCL-Hancock scheme 
shows inexact evolution of the axial component of the magnetic field B z and that this error 
converges slowly as shown in the lower right panel of figure [TT] whereas the unsplit Runge- 
Kutta scheme maintains B z = exactly to machine precision. This error is due to the failure 
of the direction split scheme to produce the an exactly divergence free representation of the 
magnetic field for the time-cen tered predictor state as discussed in the last paragraph of 



§2.41 In particular, as noted by lGardiner fc Stond (120081 ) . the evolution of the z-component 



of the magnetic induction equation reduces to dB z /dt = v z (dB x /dx + dB y /dy) so that with 
finite v z the exact evolution of B z will only be maintained for unsplit schemes with exactly 
divergence free reconstruction. 



- 49 - 




Q.5 1 1.5 2 ^ 0.2 0.4 0.6 0.8 1 



Fig. 11. — Left: Field loop advection magnetic field lines (red) over a greyscale representa- 
tion of the current density J = V x B for the initial condition (top) and after the advection 
of the loop through the periodic domain at t = 1 for both the Runge-Kutta (center) and 
MUSCL-Hancock temporal integration (bottom). Right: Evolution of the spatially aver- 
aged magnetic energy normalized to the initial analytic value B for the MUSCL-Hancock 
(top) and Runge-Kutta (center) integration. The grid resolutions shown correspond to 
iV = 32, 64 and 128 for the dash-dot, dash and solid lines respectively. Evolution of the 
spatial average of the normalized axial component of the magnetic field \B Z \ for the direc- 
tion split MUSCL-Hancock integrator (bottom). The plot shows decreasing error for the 
three resolutions N = 32, 64 and 128. 
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Table 5: MHD Loop Advection Parameters. 



p 


1 


v x 


2 


Vy 


1 


v z 


1 


p 


1 


A z 


B Q (r -R) iir<R 
otherwise 


r 


^ x - 1)2 + ( y _ 1)2 


B 


io- 3 


R 


0.3 


7 


5/3 


CFL 


0.9 (MUSCL-Hancock) 
0.45 (Runge Kutta) 
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The next set of simulations show the propagation of a cylindrical, supersonic cloud 
through a magnetized medium where the magnetic field is oriented parallel to the propagation 
of the cloud (figure O top), perpendicular to the propagation of the cloud (figure [121 middle), 
and 45° to the propagation of the cloud (figure [121 bottom). These simulations have been 
carried out using PHM spatial reconstruction, MUSCL-Hancock temporal reconstruction, 
the adapted Mar quina flux and the constrained transport evolution of the magnetic field of 



Ryu et al.l (119951 ). Each of these simulations employ constant extrapolation conditions along 



each boundary. The AMR level with the finest resolution achieves an effective resolution 
of 32 cells per cloud radius. The initial cloud density is smoothed linearly to that of the 
ambient environment and the velocity is smoothed over the outer four computational cells 
of the cloud. The figures show the result of each simulation at two evolutionary times, 
t = 2t cr in the left panel and t = 4t cr in the right panel where t cr = 1r c ^Jxjv c is the cloud 
crushing time where r c is the cloud radius, v c is the initial c loud speed and y i s the density 
contrast of the cloud against the ambient environment (see IJones et al.l (119961 ) for details). 
The density distribution is presented in gray-scale, red lines delineate the magnetic field lines 
and blue lines delineate regions of AMR-enhanced resolution. We note that the turbulent 
wake behind the 45° cloud and the associated early onset of tearing mode instability and 
magnetic reconnection pose some degree of numerical difficulty which requires the use of the 
more accurate and robust Runge-Kutta method over the faster MUSCL-Hancock method. 
The simul a tions use the same ini t ial ph ysical parameters as the simulations presented in 



Ryu et al.l (119981 ) and IJones et al.l (119961 ) which are presented in table [5l 



Because the shocked cloud simulations presented here do not utilize a moving mesh, 
the final time (t = 4thr) of the sim ulati ons is somewhat e arlier than that of the earlier 
works {t = 6t bc ) of iRyu et al.l (119951 ) and IJones et al.l (119961 ). The density distribution and 
magnetic field lines at (t = 2£& c ) show agree ment with results of th e earlier calculatio ns at 
the same evolutionary time. As pointed out in lRyu et al.l ( 1995 ) and Jones et al. ( 1996 ). the 
nonlinear evolution of the cloud depends sensitively on the exact initial perturbations which 
develop out of the geometric mismatches between the cloud and the computational grid. The 
agreement of the AMR simulations presented here with the earlier works demonstrates the 
robustness and accuracy of the divergence preserving restriction and prolongation procedures 
presented in £JU 
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Fig. 12. — Simulation results of a supersonic cylindrical cloud moving through a magne- 
tized medium with initial magnetic field oriented parallel, perpendicular and diagonal to the 
propagation of the cloud from top to bottom. The logarithm of the density distribution is 
presented in gray-scale, red lines delineate the magnetic field lines and blue lines delineate 
regions of AMR-enhanced resolution. 
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Table 6: MHD Shocked Cloud Parameters. 





ambient 


cloud 


p 


1 


10 




10 





Vy 








V z 








P 


l/ 7 


l/ 7 


P 


4 


4 


7 


5/3 


5/3 


clump radius (r c ) 


1 




grid zones per cloud radius 


32 




initial cloud position 


(2, 0) 




CFL 


0.4 
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In figure [Ti3 we re produce the light cylindrical MHD jet simulation of iRyu et al.l (119951 ) 
and iLind et al.l (119891 ) . A jet with a top-hat profile is injected into a uniformly magnetized 
environment by imposing the physical jet conditions presented in table [7] inside of r < r, 
along the z = boundary. The symmetry of the problem dictates the use of reflecting bound- 
ary conditions along the r = boundary. All other boundaries utilize constant extrapolation. 
The simulation was carried forward on an AMR hierarchy with four levels of refinement uti- 
lizing the PHM spatial integration method with Runge-Kutta temporal integratio n, the Roe 



flux u pwinding method, and the constrained transport magnetic field evolution of IRyu et al. 



Jl998h . The finest AMR level has a resolution of 32 cells per jet radius. The figures show the 
result of the simulation at 5 evolutionary times, t = 2.43, 6.57, 1 0.62, 14.76 , and 18.00. 
Note these are the same evolutionary times shown in the results of IRyu et al.l (119981 ) scaled 
to the dimensionless units defined by the parameters listed in [7| The density distribution is 
presented in gray-scale, red lines delineate the magnetic field lines and blue lines delineate 
regions of AMR-enhanced resolution. Because the simulations were carried out in cylindrical 
symmetry, the lower half of each panel has been plotted by reflecting the computational do- 
main about the symmetry axis. Barring differences in the detail of the nonlinear evolution 
of the Kelvin-Helmholtz shear flow behind the jet bow shock which arise from the exact 
nature of the numerical perturba tions present in th e simulation, the simulation results show 
excellent agreement with that of IRyu et al.l (119951 ) at each evolutionary stage shown in the 
panels of figure [131 The agreement of the AMR simulations presented here with the earlier 
works demonstrates the robustness and accuracy of the divergence preserving restriction and 
prolongation procedures presented in §3]and the adaptations thereof for cylindrical geometry 
presented in the appendix. 
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Fig. 13. — Light MHD jet simulation. The logarithm of the density distribution is presented 
in gray-scale, red lines delineate the magnetic field lines and blue lines delineate regions of 
AMR-enhanced resolution at t — 2.43, 6.57, 10.62, 14.76, and 18.00 from top to bottom. 
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Table 7: Light MHD Jet Parameters. 





ambient 


jet 


p 


0.1 


1 


V x 





20 


Vy 








Vz 








P 


l/ 7 


l/ 7 


B r 








B e {r) 





2^0.02/7 r/rj 


B z 


v/0.02/7 


v/0.02/7 


7 


5/3 


5/3 


jet radius r,j 


1 




grid zones per jet radius 


32 




CFL 


0.4 
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To confirm the robustness of the code and to demonstrate its application with radiative 
energy loss, we illustrate the propagation of a strongly radiative shock through an inhomoge- 
neous environment. The inhomogeneities in these two dimensional calculations take the form 
of a random distribution of 150 cylindrical clouds with 100 times the density of the ambient 
environment at their center. The density of the clouds follow a hyperbolic tangent smooth- 
ing function to the ambient density along the outer 50% of the radius of each cloud. The 
inhomogeneities are initially in pressure balance with the ambient environment. The upper 
and lower boundaries of the computational domain are periodic, the right boundary follows 
a constant extrapolation and the wind state is held constant in the left boundary throughout 
the simulation. The physical state of the ambient environment, the center of the clouds, and 
the ambient environment are given in tabled Figures [T4l and [151 show the results of the sim- 
ulations at several evolutionary stages (t = 0, 92, 184 and 369 yr). The simulations were 
carried forward using MUSCL spatial reconstruction, Runge-Kutta temporal i ntegration, 



the ad apted Marquina flux upwinding and the constrained transport method of iRyu et al. 



(119951 ) . An operator split energy sink source term is used to includ e the effects of radia- 



tive e nergy loss via atomic line cooling using the cooling function of iDalgarno fc McCray 



(119721 ) . The density distribution is presented in gray-scale, red lines delineate the magnetic 
field lines and blue lines delineate regions of AMR-enhanced resolution. Figure [TH shows 
the case where the magnetic field oriented parallel to the direction of wind and figure [151 
shows the case where the ambient environment is threaded with a magnetic field that is 
oriented perpendicular to the direction of an unmagnetized wind. The turbulent nature of 
the flow pattern that emerges in simulations, compounded sock compression ratios as high 
as ~ 12 which are achieved via radiative losses results in the development of persistently 
converging flow. Such flows are particularly problematic for codes that do not maintain 
the solenoidal constraint on the magn etic field as local dive rgences tend to accumulate in 
highly compressed regions of the flow (IBalsara fc Kiml 12004 ) . These simulations, therefore, 
demonstrate the robustness of the methods described in this paper for the simulation of such 
flows. 





Fig. 14. — Shock propagation through multiple clouds with magnetic field oriented parallel 
to the direction of shock propagation at evolutionary time t = 0, 92, 184 and 369 yr. The 
logarithm of the density distribution in cm~ 3 is presented in gray-scale, red lines delineate 
the magnetic field lines, blue lines delineate regions of AMR-enhanced resolution and the 
cloud diameter is used as the length unit. 
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Fig. 15. — Shock propagation through multiple clouds with magnetic field oriented perpen- 
dicular to the direction of shock propagation at evolutionary time t — 0, 92, 184 and 369 yr. 
The logarithm of the density distribution in cm" 3 is presented in gray-scale, red lines delin- 
eate the magnetic field lines, blue lines delineate regions of AMR-enhanced resolution and 
the cloud diameter is used as the length unit. 
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Table 8: Multiple Cloud Parameters. 





ambient 


cloud 


wind 


p (cm 3 ) 


100 


10 4 


100 


v x (kms -1 ) 








150 


Vy 











v z 











T(k) 


200 


0.2 


10 4 


7 


5/3 


5/3 


5/3 


cloud radius (AU) 


50 






grid zones per cloud radius 


12 






CFL 


0.4 
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5.3. Three Dimensions 



Circularly polarized Alfven waves are a n exact nonlinear so l ution to the MHD equati ons. 
We follow an approach similar to that of iGardiner &: Stond (120081 ) and iTothl (120001 ) by 
rotating the one dimensional prescription of the problem (table [9]) onto a three dimensional 
periodic grid of size v^6/2 x \/6 x \/Q with resolution iV x 2N x 2N via the rotation: 



X 




cos(a) cos(/?) 


— sin (a) — cos(a) sin(/3) 






y 




sin (a) cos(/?) 


cos(a) — sin (a) sin(/3) 




yi 


z 




sin(/3) 


cos(/3) 







(54) 



where a = arcsin and (3 = arcsin(l/v^6) so that the initial state is periodic with the 

grid and the wave-vector points along the diagonal of the computational domain. Solutions 
have been computed for one crossing time (t = 1) using monotonized-centered limited linear 
spatial reconstruction, the Roe flux function and both the unsplit Runge-Kutta and direction 
split MUSCL-Hancock temporal integrators. The left panel figure [TBI shows the convergence 
of the volume averaged norm of the L\ error vector as measured with respect to the initial 
state 

1 

e ~ 4iV 3 



\ 




\ lrj* =1 

/ j \ "i,j,k,nq 



t=0 
"i,j,k,nq \ 



(55) 



i,j,k 



for grid resolution s of N = 8 , 16,32 and 64. Consistent wit h the results of other authors on 
this test problem (lT6thll2000l ; IGardiner fc Stondl2005l . 120081 ) . the solution error arises mainly 
from the magnetic field components that are transverse to the wave propagation direction. 
The right panel of figure [16] shows the steady convergence of the amplitude of the transverse 
components of the magnetic Bt = 
space 



B^ + B 2 zX field after rotating back into the unrotated 



Xx 




cos(a) cos(/3) 


sin(a) cos(/3) sin(/?) 




X 


y\ 




— sin(a) 


cos(a) 




y 






— cos(a) sin(/3) 


— sin(a) sin(/3) cos(/?) 




z 



(56) 



for the unsplit Runge-Kutta case. Both integration techniques converge in a manner con- 
sistent with the expected second order accuracy with least squared power-law induces e oc 
jy-1.92 f or MUSCL-Hancock scheme and e oc N~ 2A9 for the Runge-Kutta scheme. 
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Fig. 16. — Circularly polarized Alfven wave in three dimensions: Left: convergence of 
the L\ error for the direction split MUSCL-Hancock temporal reconstruction (dotted) 
and the unsplit Runge Kutta temporal reconstruction (solid), Right: transverse compo- 
nent of the magnetic field (right) after one grid crossing (t = 1) for grid resolutions 
N = 8 (dot) , 16 (dash-dot , 32 (dash) and 64 (solid) . 
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Table 9: Circularly Polarized Alfven Wave Parameters. 



p 


1 


Vl 





V2 


a sin(27ra;i) 


v 3 


a cos(27ra;i) 


P 


1 


B 1 


1 


B 2 


a sin(27ra:i) 


Bz 


a cos(27rxi) 


a 


0.1 


R 


0.3 


7 


5/3 


At 


l/(3iV) (MUSCL-Hancock) 




l/(6iV) (Runge Kutta) 
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6. Conclusion 

The staggered grid constrained transport schemes described in this paper enable the 
application of high resolution shock capturing methods to magnetized flow. In this paper 
we have demonstrated that a wide cross section of high resolution shock capturing schemes 
for general conservation laws may be adapted for magnetized flow while preserving the 
divergence free constraint on the magnetic field topology exactly by conserving the surface 
integral of magnetic flux over each computational cell in an upwind fashion. The use of 
such schemes on multi-resolution AMR grids is encumbered by the requirement that the 
prolongation and restriction steps preserve the divergence free topology of the magnetic fields. 
In this paper we have described the application of prolongation and restriction operators 
which maintain such topologies to machine precision. 

The numerical schemes discussed here have been implemented and tested in the As- 
troBEAR adaptive mesh refinement code. The code utilizes a modular design, enabling 
the user to choose from various methodologies to tailor the numerical integration strategy 
to the requirements of the application at hand. The robustness of this approach to high 
resolution, shock capturing MHD on AMR grid structures, and relative advantages of the 
various numerical schemes implemented in the code are demonstrated in the context of sev- 
eral numerical example problems. The description of the numerical schemes presented in this 
paper provides a concise recipe for their implementation which will enable the reproduction 
of these outcomes by other researchers and the interpretation of future works derived from 
the AstroBEAR code. 
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A. Cylindrical Axisymmetry 

The conservative update procedure for integrating the ideal MHD equations $2] may be 
readily extended to the case cylindrical axisymmetric flows via a change in the coordinate 
variables and the addition of a source term to equation [2j In particular, the coordinate 
subscripts transform as as (x, y,z),—* (r, 9, z) and the source term becomes: 

S = -~ [pv r , pvl - + B e - pvg, 2(pv r v e - B e B r ), pv r v z - B z B r , 

v r (E + P + B 2 /2) - B r (B • v), 0, 0, v r B z - v z B r f . (Al) 

Due to axisymmetry, all differential terms in 9 vanish. We handle the geometric source term, 
S separately from the conservative update of the homogeneous part of the system using 
an operator split approach. The source term step ^ = S is integrated via a fourth-order 



Rose nbrock integration scheme for stiff systems of ordinary differential equations (IPress et al. 



19921 ) using an adaptive time-step to maintain the accuracy of the solution to a user specified 
level (usually 1 part in 10 4 ). Symmetry dictates that such simulations be carried out in a 
half-meridional plane (r > 0) and that reflecting boundary conditions be applied along the 
r = plane. 

The constrained transport of §2.41 and prolongation / restriction of §H of grid-interface 
magnetic field components, B r and B z may also be readily adapted for cylindrical axisym- 
metric flow. Re aders interested in the extention to more complicated geometries may refer 



to the work of iBalsaral (120041 ) which derives divegence free reconstruction procedures in 
several three dimensional curvilinear geometries and on tetrahedral meshes. In cylindrical 
axisymmetry, the solenoidality constraint on the magnetic field takes the form: 

„ „ ldrB r dB z /A . 

V-B = -— ^ + — ^ = 0, (A2) 
r or oz 

and the component of the curl of the magnetic field orthogonal to the symmetry plane takes 
the form: 

T „, dB r dB 7 , . „. 

i v x B '» = ar - (A3) 

These operators take the same form as the Cartesian case under the change of variable 
B r — > rB r . Therefore CT update formulae of §2.41 and the prolongation and restriction 
forulae of §@]are written for the case of axisymmetry by replacing (f x , f z ) — > (f r , f z ), E y — > 
Eg, (B X ,B Z ) — > (rB r ,B z ), and (b x ,b z ) — > (rb r ,b z ). The CT integration procedure (equations 
I3TT) can be written for axisymmetric geometry by replacing B z — > rB z and E y — > tEq. 
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B. Pseudo-code listing for the AMR engine. 

SUBROUTINE AMR(level , dt) 
IF(level=0) 

nsteps = 1 

Set Ghost (level) 
ELSE 

nsteps = r 
END IF 

DO n = 1, nsteps 
Distribute (level) 
IF (level < MaxLevel) 

Grid Adapt (level + 1) 

Set Ghost (level+1) 
END IF 

Integrate (level , n) 

IF(level< MaxLevel) AMRQevel + 1, dt/r) 
END DO 

IF(level > 1) Syncronize to Parent (level) 
END SUBROUTINE AMR 
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